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Abstract. We study the evolution of circumbinary disks surrounding classical T Tau stars. High resolution numer- 
ical simulations are employed to model a system consisting of a central eccentric binary star within an accretion 
disk. The disk is assumed to be infinitesimally thin, however a detailed energy balance including viscous heating 
and radiative cooling is applied. A novel numerical approach using a parallelized Dual-Grid technique on two 
different coordinate systems has been implemented. 

Physical parameters of the setup are chosen to model the close systems of DQ Tau and AK Sco, as well as 
the wider systems of GG Tau and UY Aur. Our main findings are for the tight binaries a substantial flow of 
material through the disk gap which is accreted onto the central stars in a phase dependent process. We are 
able to constrain the parameters of the systems by matching both accretion rates and derived spectral energy 
distributions to observational data where available. 
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1. Introduction 

Observations of main sequence stars have shown that 
about 60% belong to binary or multiple systems 



( Duqucnnoy fc Mayor 1991 ). However, recent observations 
of pre-main sequence stars seem to indicate that nearly all 



stars are b orn in multiple systems (R eviews in Mathicu 



ct al. 2000; Zinncckcr fc Mathieu 2001) 



Apart from many spectroscopic binary systems with 
surrounding disks like DQ Tau and AK Sco, there have 
been only two observations of directly imaged wide opti- 
cal binary systems with circumbinary disks, namely GG 
Tau dDutrcy ct al. 1994| ) and UY Aur ( |Duvcrt et al. 1998j ). 
The observations of the spectroscopic binaries usually in- 
clude emission data from the infrared and optical part 
of the spectrum. Observations of the luminosities of the 
boundary layer are typically used to estimate accretion 



ct al. 1W5J 



rates o rto the binary stars ^Gullbrmg ct al. lyysfHartigan ^ L u t, ow 1994) 



In a system consisting of a stellar binary surrounded 
by an accretion disk there is a continuous exchange of an- 
gular momentum between the binary and the disk. As the 
binary is orbiting with a shorter period than the disk ma- 
terial it is generating positive angular momentum waves 
in the disk which, upon dissipating, deposit their energy 
and angular momentum in the disk. Hence, there is a con- 
tinuous transfer of angular momentum from the binary to 
the outer disk. Upon receiving this angular momentum, 
the inner disk speeds up, and the material is receding from 
the binary, leaving a region of very low density around the 
central binary, a so called inner gap. At the same time, the 
loss of angular momentum leads to a secular shrinkage of 
the semi-major axis of the binary. The width and form of 
the gap depend on disk parameters such as temperature 
and viscosity and on the orbital parameters (eccentricity) 
of the binary as well (Artymowicz et al. 1991; Artymowicz 



Bate & Bonncll (1997) model the early phase of a 
molecular cloud collapsing onto a protobinary system with 
an SPH code and derive properties of the created cir- 
cumbinary and circumstellar disks. We are interested in 
the following phase of the evolution when the envelope 
has been cleared. By modeling the evolution we compare 
fully developed circumbinary and circumstellar disks and 
their properties with observational data. 



First calculations dealing with the evolution of cir- 
cumbinary disks so far have been using the numerical 
method of smoothed particle hydrodynamics. Applying 
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such a scheme Artymowicz & Lubow (1996) could demon- 
strate that for sufficiently high viscosities and tempera- 
tures material from the disk is able to penetrate into the 
inner gap region of the disk to become eventually accreted 
by one of the binary stars. This process which is facilitated 
by a gas flow along the saddle points of the potential leads 
to a pulsed accretion whose magnitude and phase depen- 
dence for eccentric binaries was estimated. However, when 
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using only a small number of particles all having the same 
mass only a limited resolution is achievable. 

Later Rozyczka & Laughlin (1997) used a finite dif- 
ference method in attacking this problem. Using a con- 
stant kinematic viscosity coefficient and eccentric binaries 



they also found, in agreement with Artymowicz & Lubow 
(1996|)] a pulsed accretion flow across the gap, with most 
of the accretion occurring onto the lower mass secondary 
star. 

In this paper we extend these computations and solve 
explicitly a time dependent energy equation which in- 
cludes the effects of viscous heating and radiative cooling. 
We aim at the structure and dynamics of the disk as well 
as the gas flow in the close vicinity of the binary star. To 
that purpose we utilize a newly developed method which 
enables us to cover the whole spatial domain. For the first 
time we perform long-time integration of the complete sys- 
tem covering several hundred orbital periods of the binary 
and compare the properties of the evolved systems with 
observational data such as spectral energy distributions 
in the infrared and optical bands and accretion rates esti- 
mated from luminosities. 

In the next section we layout the physical model fol- 
lowed by some remarks concerning important numerical 
issues (Sect. |J). We describe the setup of the various nu- 
merical models in Sect. ^. The main results are presented 
in Sect. ^ and our conclusions are given in Sect. |^. 

2. Physical Model 

2.1. General Layout 

The system consists of a stellar binary system and a sur- 
rounding circumbinary accretion disk. We assume a ge- 
ometrically thin disk and describe the disk structure by 
means of a two-dimensional, infinitesimal thin model in 
the z — plane with the origin at the center of mass 
(COM) of the binary. We use vertically averaged quanti- 
ties, such as the surface mass density 



£ = 



pdz, 



where p is the regular three-dimensional density. We work 
with two different coordinate systems. Firstly, a cylindri- 
cal coordinate system (r, ip) centered on the COM of the 
binary stars. As such a coordinate system does not al- 
low for a full coverage of the plane, we overlay the center 
additionally with a Cartesian (x, y) grid (see Sec. [T^j- 

The gas in the disk is non-self-gravitating and is or- 
biting a binary system with stellar masses Mi and M2 
which are modeled by the gravitational potential of two 
solid spheres. The binary is fixed on a Keplerian orbit 
with given semimajor axis a and eccentricity e. Hence, we 
neglect self-gravity and the gravitational backreaction of 
the disk onto the stars. This is justified because the total 
mass of the disk Md within the simulated region, which 
extends up to several tens of semi-major axis a of the bi- 
nary, is typically much smaller than that of the central 



binary. Hence, the timescale for a change in a and e is 
indeed much longer than the timescales considered here 
(Artymowicz ct al. 1991). For the majority of the compu- 
tations we employ an inertial coordinate system, however 
test calculations were performed in a system corotating 
with the binary as well. 



2.2. Equations 

The evolution of the disk is given by the two-dimensional 
(x, y or r, tp) evolutionary equations for the density S, the 
velocity field u = (u x , u y ) = (u r , it v ), and the temperature 
T. In a coordinate-free representation the equations read 



AT 

— +V.(E«) = 



dEu 

9£etot 
dt 



(1) 
(2) 



V • (Su • u) = -VP - SV$ + V • T 
+ V-[(Ze tot + P)u] = -Eu4-V-(f-ii-T)(3) 



Here e tot = l/2u 2 + e t h is the sum of the specific in- 
ternal (eth = c v T) and kinetic energy. P is the vertically 
integrated (two-dimensional) pressure 



P = 



RET 



with the midplane temperature T and the mean molecular 
weight p. which is obtained by solving the Saha rate equa- 
tions for Hydrogen dissociation and ionization and Helium 
ionization. 

The gravitational potential <f>, generated by the binary 
stars, is given by 

GAL 



E 

i=l,2 



GMi(ZRl-\r-n\ 2 ) 
2Rl 



for li- 



fer Ir 



> R 



(4) 



< R 



where r, are the radius vectors to the two stars, and R is 
the stellar radius assumed to be identical for both stars. 
The second case in Eq. (^) gives the potential inside a star 
with radius R having a homogeneous (constant) density. 

2.3. Viscosity and Disk Height 

The effects of viscosity are contained in the viscous stress 
tensor T. Here we assume that the accretion disk may be 
described as a viscous medium driven by some internal 
turbulence which we approximate with a Reynolds ansatz 
for the stress tensor. The components of T in different 
coordinate systems are spelled out explicitly for example 
in Tassoul (1978). A useful form for disk calculations con- 



sidering angular momentum conservation is given in Kley 
(19991 ) • 

For the kinematic shear viscosity v we assume in this 



work an a- model (Shakura & Sunyaev 1973) in the form 



aCsH 



(5) 
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(6) 



where c s and H are the local sound-speed and vertical 
height, respectively. The local vertical height H(r) is com- 
puted from the vertical hydrostatic equilibrium 

9$ 1 dp 

dz p dz 

where p and p are the standard three dimensional density 
and pressure. Substituting Eq. ([|) for the potential with 
|r > i?*, assuming an ideal gas law p = RpT/p,, and 
integrating over z yields 

p = p e~(^) (7) 

where po = S/[(27r) 1 ' 2 H] is the midplane density, and H 
is the vertical height given by 



H(r) = 



GAL 



= 1,2 C s \ r - r * 



i=1.2 



H-\r) (8) 



which can be split into the single star disk heights given 
by Hi(r) = c s ^J ^jj^-- The sound speed is given here by 
c s = RT / p. We assume a vanishing physical bulk viscosity 
C, but consider it in the artificial viscosity coefficient (see 
Kley 1999|). 



2.4. Radiative Balance 

The influence of radiative and viscous effects on the disk 
temperature are treated in a local balance. The balance 
equation for the internal and radiative energy considering 
only these two effects reads 



Q 



diss 



; racl 



(9) 



d (Ec v T + aT 4 ) 
dt 

where Qdiss and Qrad denote the viscous dissipation and 
radiative losses, respectively. For the viscous losses we use 
the vertically averaged expression 

(10) 



u-VT 



-Tr (T 2 ) 



; rad 



(11) 



2zZv 

For the radiative transport 
dF z 
dz 

where Fo is the flux vector in the z = plane. Here 
we consider only the losses in the vertical direction, i.e. 
Fq = 0, a standard approximation in accretion disk the- 
ory. Integration over the vertical direction yields 

Qrad = 2F ra d = 2<7-qT^ s (12) 

with the local effective (surface) temperature T e s which is 
related to the midplane temperature T through 

-1/4 



'off 



1/4 T 



2 4 



T 



(13) 



using an interpolated opacity n(po, T) adapted from Lin & 
Papaloizou (198E). Note that this formulation of radiative 
transfer is a very good approximation only in the optically 
thick regime. 



3. Numerical Issues 

In order to study the disk evolution, we utilize a finite 
difference method to solve the hydrodynamic equations 
outlined in the previous section. As we intend to model 
possible accretion flow onto the binary stars and to achieve 
a very good resolution in their vicinity, we use a Dual- 
Grirf-technique . 

3.1. General 

The code is based on the hydrodynamics code RH2D 
suited to study general two-dimensional systems includ- 
ing radiative transport (Kley 198£). RH2D uses a spa- 
tially second-order accurate, mixed explicit and implicit 
method. Due to an operator-splitting method it is semi- 
second order accurate in time. The advection is computed 
by means of the second order monotonic transport al- 



gorithm, introduced by van Leer (1977), which guaran- 
tees global conservation of mass and angular momentum. 
Advection and forces are solved explicitly, and the viscos- 
ity and radiation are treated implicitly. The formalism for 
application to thin disks in (r — ip) g eometry has been de- 
scribed in detail in ( Kley 1998 , 1999 ). Here thi s sch eme is 
extended to allow for a Dual-Grid system (see ^2) and a 
more detailed energy balance Eq. (||). 

As the calculations have shown, the system behaves ex- 
tremely dynamically with very strong, moving gradients 
at the inner edge of the disk. To model such a situation, 
the viscosity has to be treated implicitly to ensure numeri- 
cal stability. The coupling of u r and u v requires a solution 



method similar to that in Kley (1989). For stability pur- 
poses an additional implicit artificial viscosity has to be 
included in the computations ( Kley 1999] ). 

This code has been employed successfully in related 



simulations dealing with embedded planets in disks (Kley 
|1999| , |2000| ). The numerical details of the present simula- 
tions are presented in Giinther (2001). 



3.2. Dual-Grid technique 

The Dual- Grid technique employed here combines two de- 
sirable features of thin disk computations: a) An ideally 
suited coordinate system (r — tp) in the outer bulk parts 
of the disk which ensures conservation of angular momen- 
tum, and b) Coverage of the region near the central binary 
star allowing the determination of the mass flow onto the 
stars. 

In the main part the domain is covered by a cylindrical 
grid. On this (outer) grid the equations arc formulated in 
polar coordinates, and the solution technique guarantees 
exact conservation of global angular momentum. The con- 
servation property is absolutely crucial in obtaining phys- 
ically reliable results for long term disk evolution ( Kley 
|1998D . 

However, such a grid cannot be extended to the very 
center and leaves a hole in the middle. For typical com- 
putations this does not pose serious problems, but in this 
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Fig. 1. Grid system of a typical calculation. Shown is the 
grid structure of the Cartesian and cylindrical grid at the 
center of the computational domain. Circumstellar disks 
of the system can be seen at the top and bottom edges of 
the picture. 



case we are dealing with eccentric binary systems in the 
middle. Hence, the stars or their circumstellar disks may 
easily cross the inner grid boundary causing strong irreg- 
ularities. In any case a central hole does not allow for 
an accurate determination of the individual mass accre- 
tion rates onto or a possible overflow between the stars. 
Therefore, using a Cartesian grid centered on the origin 
enables an accurate calculation of the mass flows in the 
vicinity of the stars. Applying this technique, a high lo- 
cal resolution can be obtained in the center as well. An 
example of such a grid structure is displayed in Fig. [I]. 

Related numerical schemes using a nested grid tech- 
nique, though not in different coordinate systems, have 
been applied in astrophysical simulations by a number of 



authors (eg. [Ruffert 1992| ; |Yorkc et al. 1993[ ). 

The necessary equations (in Cartesian and cylindri- 
cal coordinates) are then integrated, independently and 
in parallel, on the two grids, and after each timestep the 
necessary information in the overlap region has to be ex- 
changed. Restrictions are imposed on the time step for 
stability reasons, the Courant-Friedrichs-Lewy (CFL) con- 
dition must be fulfilled during each integration, on each 
grid. 

Test calculations using a single Cartesian grid covering 
the whole domain have shown that the non-conservation of 
angular momentum yields serious deviations from known 
analytical results in particular for longer integration times. 

As the grids possess a different geometry, exact con- 
servation of momenta and energy cannot be guaranteed 
in the conversion/interpolation process which uses linear 
interpolation to exchange density, temperature and veloc- 
ity. Exchanging density, internal energy and momentum 
was also implemented and leads to more correct conserva- 
tion of momenta and energy. However, the advantage of 




Fig. 2. Radial shock in first (infalling) phase right af- 
ter crossing the grid boundaries. Surface density is color 
coded, the grid structure is shown. 

a complete coverage of the physical space in between the 
two stars clearly outweighs this disadvantage. 

3.3. Tests 

We ran several tests to check the accuracy and numerical 
stability of the grid overlap region, including an infalling 
shock and a Sedov like explosion. For these tests no viscos- 
ity nor radiative effects are applied. All tests show that the 
grid overlap technique is sufficient for handling problems 
with radial symmetric features such as accretion disks. 

The shock experiment initial setup is a step function 
of density and temperature with the floor starting at ro 
(about 10 grid cells outside of the inner boundary of the 
cylindrical grid) and covering the whole Cartesian grid: 



£(r) = S n 



S 6(r 



ro) 



(14) 



In Fig. g you can see the grid setup and the shock in the 
initial infalling phase just after crossing the grid bound- 
aries. The shock front remains radially symmetric and no 
reflections are visible. After the shock compresses into a 
single point it switches to the outflowing phase. You can 
see the shock front propagating outwards just after cross- 
ing the grid boundaries in Fig. ^. Apart from unavoidable 
grid effects near the center the shock front remains nicely 
circular and again neither disturbances nor reflections are 
visible at the grid interface. 

4. General model design 

The main goal of this study is the investigation of the 
characteristic features of the hydrodynamic flow of the 
circumbinary disk. 
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Fig. 3. Radial shock in second (explosion) phase right af- 
ter crossing the grid boundaries again. Surface density is 
color coded, the imaged domain is the same as in Fig. 0. 



From now on, we refer to non-dimensional units. All 
the lengths are expressed in units of the semi-major axis 
a of the binary. This is constant because the time scale 
on which the binary orbital elements change due to mass 
and angular momentum transfer is much longer than the 



dynamical time scales considered here (see Artymowicz & 



Lubow 2001). Masses are in units of the solar mass and 
time is given in units of the binary's orbital period. 

The whole azimuthal range of the disk is taken into ac- 
count by considering a computational domain represented 
by 2tt x [0, r ut]) where typically r out = 10 in units of 
the binary separation. The cylindrical grid covers a re- 
gion 2ir x [rid, r ou t] with r- 1B = 0.1 represented by up to 
256 x 256 meshpoints, and the central Cartesian grid has 
up to 128 x 128 meshpoints covering 2 r; n x 2 r ln with an 
additional small overlap. 

4.1. Accretion onto the Stars 

As will be seen, matter flows from the disk through the 
gap and circulates the stars in small circumstellar disks. 
From those it can be accreted onto the stars. This stellar 
accretion is accounted for by removing some mass from 
a region close to the stars. This removal is performed on 
both grids dependent on the star positions. The removed 
mass is not added to the dynamical mass of the stars, but 
is just monitored. 

The region R out of which matter is removed is given as 
a fraction (0.2) of the Roche lobe size for the wide systems 
and as the stellar cores for the tight spectroscopic systems. 
The accretion process typically involves only a few grid 
cells on each side of the stars for the wide systems and 



the whole resolved stellar core (up to 11x11 grid-cells) for 
the tight systems, making it a locally confined process. 

Two methods for deciding whether a grid cell is tar- 
geted for accreting a part of its mass were implemented. 
Accreting a constant fraction / per time of the mass ex- 
ceeding a defined minimum density S m ; n leads to a very 
smooth accretion process. The accreted mass per iteration 
is computed by 

M acc = min(1.0,/At)max(0.0,Sy-S mi „)AVol ?J (15) 

ijeR 

Limiting the density in the accretion region R to an ar- 
bitrary value S max results in a more realistic structure of 
the circumstellar disks and to more correct spectral energy 
distributions of the circumstellar disks. 



(16) 



M acc = max(0.0, E« - S max ) A Vol, 

i,j€R 

This accretion process is not smooth, but consists of ac- 
cretion events because exceeding the maximum density is 
a discrete process and not distributed evenly over time. 
Nevertheless the time-averaged accretion rates are the 
same for both methods. 

The inferred accretion rates in any event compare fa- 
vorably with the average disk accretion rate, as estimated 
from the radial mass fluxes through different surfaces, and 
the evolution of the total mass found in the circumstellar 
disks. 

4.2. Spectra generation 

Spectral energy distributions are generated decoupled 
from the dynamic evolution of th e disk as a post- 
processing step. Here the model from Adams ct al. (1988 ) 
is extended to work on data generated by our numerical 
simulation. The flux at the observer is obtained by inte- 
grating over the disk surface assuming local black-body 
radiation 



F v = 



COS I 



B v [T(r,<p)] 1-e 



,-t{t,<p) 



rdip dr(17) 



where i is the inclination of the disk, D the distance to 
the observer and r the line-of-sight optical depth through 
the disk, r can be obtained using the opacity k, t(t, ip) = 
«s(r,y) rp j g ^ c sur f acc ^jgjj temperature which has to 

COS I ^ 

be estimated using the optical depth of the disk obtained 
from an opacity table. 

The integration is limited to regions outside of the stel- 
lar cores as temperatures and densities from inside the 
cores can be nothing but wrong. The star emission is ac- 
counted for by adding a black-body spectrum from an 
appropriate effective star temperature. 



.-d 



4.3. Initial and boundary conditions 

The initial density distribution is proportional to 
where d is either set to 0.5 which is the equilibrium den- 
sity distribution for an accretion disk around a central star 
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Fig. 4. Initial gap profile composed out of Y.qt 1 and the 



gap function Eq. 
with Sq = 1. 



jlq). The figure shows a gap at 



= f 



having a constant kinematic viscosity, or given by a fit of 
the generated spectral energy distribution to the one ob- 
served. However, we superimpose to this an axisymmetric 
gap around the COM of the two stars, given by an ap- 



proximate gap radius r- gap ( Artymowicz fc Lubow 1994 ) 
and a gap function 



gap •- 



1 



(18) 



Figure [| shows a typical surface density at t = 0. The 
initial velocity field is that of a Keplerian disk orbiting 
the center of mass of the two stars. 

The gravitational potential in Eq. (Q) is phased-in 
smoothly from an initial central potential for the total 
mass, M\ + M2, during an initialization phase lasting typ- 
ically ten orbital periods. This allows the system to adjust 
smoothly from the semi-stationary setup of a single cen- 
tral star to the binary star situation. 

For the boundary conditions, periodicity is imposed at 
ip = and f = 2 tt. We set outflowing boundary con- 
ditions at the outer radial border (r out ), i.e. the radial 
velocity gradient at -R max is zero, and truncated to posi- 
tive values. Angular velocity there equals the unperturbed 
initial Keplerian value. 

For single grid calculations with a cylindrical grid we 
use zero gradient boundary conditions for both velocity 
components and the density at the inner border. 

5. Results 

Before discussing the results concerning particular astro- 
nomical objects, we would like to present the details of a 
specific test case of an equal mass binary star on a circu- 
lar orbit. Such a system is well suited to check important 
properties (symmetries, stability) of the numerical method 
used. Also general properties of accretion disks in binary 
systems can be derived from such a test case. 




Fig. 5. Evolution of an equal mass binary after 90 or- 
bital periods. Color encoding is logarithmic surface den- 
sity. Both the inner part of the circumbinary disk and the 
circumstcllar disks are visible. 



5.1. Equal Mass Binary 

The equal mass binary on a circular orbit serves as a test- 
ing platform for numerical stability and accuracy. In a co- 
ordinate frame corotating with the binary one expects af- 
ter some transient time a quasi-stationary situation, which 
is symmetric with respect to the line connecting the stars. 
As the viscosity is solved implicitly by solving a large 
matrix system iteratively, exact symmetry cannot be pre- 
served numerically. Also usual density contrasts are of the 
order of 10 3 for numerical floor to circumstellar disks and 
from circumstellar disks to the circumbinary disk. Hence, 
this test serves as an excellent example for estimating the 
ability to preserve symmetry over long time scales. 

In Fig. ^|you can see a circumbinary disk after 90 orbits 
evolution time in almost perfect quasi-stationary state. 
The deviation from line-symmetry is small. Only with very 
high resolution simulations we observe spiral features in 
the outer regions of the disk, while spiral arms originating 
from the inner gap of the circumbinary disk flowing onto 
the circumstellar disks can be observed even for low reso- 
lutions. Those spiral arms feeding the circumstellar disks 
remain stable for the whole period of a non-eccentric bi- 
nary. 

The accretion rates for both accretion mechanisms 
onto the two stars are nearly identical and do not show any 
significant variations for different numerical parameters. 
Thus, our implementation of the accretion mechanism is 
robust. 

Changing model parameters such as So, gap width, 
or disk size we obtain characteristic variations of the re- 
sulting spectral energy distribution of an evolved disk. 
Figures ^|-^ show dependencies of these parameters on the 
spectral energy distribution that are used to fit parameters 
of the spectroscopic systems to match observed spectra. 
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no initial gap 
1 AU 
2.5 AU 
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10 AU 




1e+15 





AK Scorpii 


DQ Tauri 


p 


13 d . 60933 ±0 d . 00040 


15 d . 8043 ±0 d . 0024 


e 


0.469 ± 0.004 


0.556 ±0.018 


q 


0.987 ± 0.007 


0.97 ±0.15 


i 


63° 


23° 


a 


0.16 AU 


0.135 AU 


T* 


6500 K 


4000 K 


i?* 


2.2 R 


1.785 R Q 


M* 


1.5 M Q 


1.3M q 


d 


152 pc 


140 pc 


M d 


0.005 M 


0.001 M 


R d 


40 AU 


50 AU 



Table 1. Parameters for AK Sco and DQ Tau. 



Fig. 6. Dependency of the gap width on the SED. 
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Fig. 7. Dependency of the disk size on the SED. 
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Fig. 8. Dependency of E on the SED. 



5.2. Spectroscopic Systems 

Hereafter we discuss models for the two spectroscopic 
young binary stars DQ Tau and AK Sco which have a 
separation of only a fraction of an AU. The physical pa- 



rameter for the systems are taken from Andersen et al 
(19890 [Favata et al. (199% |Jensen fc Mathicu (1997| ) for 
AK Sco and [Mathieu et al. (1997|) for DQ Tau. The rele- 




Fig. 9. AK Sco circumbinary disk after 82.0 orbital pe- 
riods in apastron. Color coding is log(S), the size of the 
stars reflects the actual stellar radii, the length scales are 
in AU 

vant information for the simulations has been summarized 
in Table [j], disk masses Md are integrated densities over 
the computational domain specified by the disk radii R^. 
Images showing the disk configuration in apastron and 
periastron phase are shown in Figs. ^| and |^| 

Generally, for eccentric systems we observe spiral fea- 
tures in the outer circumbinary disk only for very high 
resolution simulations, while again spiral arms feed the 
circumstcllar disks even with low resolutions. Opposed to 
the non-eccentric systems the spiral arms develop during 
apastron phase and are torn off together with the outer 
parts of the circumstellar disks after periastron phase for 
high eccentric systems. This tearing off is caused by raising 
centrifugal forces acting on the circumstellar material. 

For both systems we compared the accretion rates re- 
sulting from the simulation with accretion rates derived 



from of observations as done by Hartigan et al. (1995) and 
pullbring et al. (1998| ) (see Table |2j). As no data was avail- 
able for AK Sco, a comparison with observational data was 
not possible, but as both spectroscopic systems have sim- 
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Fig. 10. DQ Tau circumbinary disk after 85.5 orbital pe- 
riods in periastron. Color coding is log(S), the size of the 
stars reflects the actual stellar radii, the length scales are 
in AU. 





Hartigan 


Gullbring 


simulation 


DQ Tau 


0.50 ■ 10"' 


0.60 ■ 10" M 


0.62 • 10" 8 


AK Sco 


n.a. 


n.a. 


0.83- 10~ s 


GG Tau 


0.20 ■ 10~ b 


0.175 ■ 10"' 


0.11 ■ 10" 8 


UY Aur 


0.25 • 10~ B 


0.656 ■ 10"' 


0.5 • 10"' 



Table 2. Accretion rates in M Q yr~ 1 derived from obser 
vational data from Hartigan ct al. (1995| ) and Gullbring 
et al. (199S ) compared to rates as resulting from our sim- 
ulations. Our rates are for the primary component of the 
binaries, for the other data we assume it is an averaged 
accretion rate for both stars. 



ilar system parameters, accretion rates of the same order 
of magnitude can be expected. 

As can be seen from Figs, [ll] and [l2] the accretion 
process happens synchronized with periastron phase of 
the binaries. This is in agreement with numerical simu- 
lations from Artymowicz fc Lubow (1996 ) and Rozyczka 
& Lau jdilin (1997 ) and with obs ervations which show ex- 
cess luminosity at these phases ( Mathicu ct al. 1997 ). 

From our simulations we can derive an averaged ac- 
cretion rate of 0.83 • 10~ 8 Moyr" 1 for the primary, 0.82 • 
lO~ 8 M yr" 1 for the secondary of AK Sco and 0.62 • 
lO" 8 M yr" 1 for the primary, 0.51 ■ lO" 8 M yr~ 1 for the 
secondary of DQ Tau. These accretion rates are in good 
agreement with the observational rates. 

It is expected that the luminosity increases in pe- 
riastron because tidal compression of the disk leads to 
a temperature rise with subsequent luminosity increase. 
We tried to generate luminosity curves to show that in- 
deed excess luminosity is produced at periastron phase, 
but unfortunately those curves show both brightening 
and darkening dependent on the implemented accretion 
mechanism. At the very low averaged accretion rates of 
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Fig. 11. Time-dependent accretion rate of the AK Sco 
system. The orbital phase curve schematically shows the 
distance of the two stars. 
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Fig. 12. Time-dependent accretion rate of the DQ Tau 
system. The orbital phase curve schematically shows the 
distance of the two stars. 



M ~ 7 • lO~ 9 M0yr _1 the disk is very optically thin. In 
that regime the implemented radiative loss mechanism is 
not very accurate and thus these results are not meaning- 
ful. 

For the spectroscopic binaries we can generate spec- 
tral energy distributions out of the simulation data which 
can be made to match the observed ones by manually ad- 
justing initial model parameters such as disk mass, gap 
width and density profile d. For DQ Tau Fig. [l3] shows 
the spectral energy distribution composed out of the disk 
and the star, Fig. [li] shows the one for the AK Sco sys- 
tem. The effective stellar temperatures used in the model 
have been adjusted to match the observational data. They 
differ from the quoted values in the literature in Tabic [l]. 
This may be due to difficulties in separating individual 
stellar and disk contributions to the total flux in models 
and observations. Also occultation is not accounted for on 
adding the star spectrum to the disk. 
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190 AU 


67 AU 
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4000 K 


3800 K 
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2.6 R Q 
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0.95 M Q 


0.65 M 
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150 pc 


M d 


1.1 M Q 


0.15 M Q 


R& 


2100 AU 
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Table 3. Parameters for UY Aur and GG Tau. 



Fig. 13. Spectral energy distribution of DQ Tau (d = 2.0) 
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Fig. 14. Spectral energy distribution of AK Sco (d = 1.5 



crosses show observational data taken from Jensen & 



Mathieu (1997Q 



5.3. Wider Systems 

There are only two well known systems where the binary 



and th e dis k have been imaged dire ctly, GG Tau ( Dutrey 



994 



L998 



Guillotcau ct al. 1999|) and UY Aur flClose 



Duvcrt et al. 1998| ). These systems have a 



et al. 
et al. 

much wider separation of tens to over 100 AU, a mass 
ratio different from unity and show consequently a quite 
different behavior. In Table [| the system parameters used 
for the simulations are shown. 

From our simulations we derive averaged accretion 
rates for both UY Aur and GG Tau that are compa- 
rable with the observed quantities in Table |2| In the 
case of UY Aur we get 0.5 • 10 _7 M Q yr _1 for the pri- 
mary and 0.25 • 10 _6 M Q yr _1 for the secondary, respec- 
tively. Simulations of GG Tau lead to accretion rates of 
0.11 • 10" 8 Mgyr" 1 for the primary and 0.33- 10~ 9 Moyr" 1 
for the secondary. As expected from the difference of the 
disk masses (0.15 Mq vs. I.IMq) the accretion rate for 




Fig. 15. UY Aur circumbinary disk after 20.5 orbital pe- 
riods. Color coding is log(S), the length scales are in AU. 



GG Tau is one order of magnitude less than the one for 
UY Aur. 

Accretion is reinforced at periastron phases as for the 
spectroscopic systems, but events are an order of magni- 
tude lower for GG Tau and nearly vanish for UY Aur due 
to the lower eccentricity of the systems. Also both systems 
show more phase dependent accretion for the secondary 
than for the primary. 

In Figs. 15 and |l^ you can see the inner parts of the 
simulated circumbinary disks. The circumstellar disks re- 
main intact over the whole period. It was not possible to 
generate meaningful spectral energy distributions of the 
disks as they are too cold for these wide systems. Also 
observational data is lacking for these systems. 



6. Conclusion 

With the present calculations we have modeled the accre- 
tion flow onto the stars within a circumbinary disk in more 
detail using more realistic physics, much higher resolution 
and a much longer time integration. Our main achieve- 
ments may be summarized as follows: 
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Fig. 16. GG Tau circumbinary disk after 56.4 orbital pe- 
riods. Color coding is log(S), the length scales are in AU. 



1. Development of a new Dual-Grid technique combin- 
ing two coordinate systems, which is ideally suited for 
modeling planar accretion disks including their central 
objects. 

2. Performing long-time integrations of binary stars and 
circumbinary disks, including a detailed vertical energy 
balance. 

3. For close binaries with separations of only a fraction of 
an AU we find a relatively large accretion rate crossing 
the gap of the order of 10 _8 MQyr^ 1 . 

4. For wide binaries (a — 10 — 200 AU) the accretion rate 
is comparable to those of the close binaries only be- 
cause of the large disk masses of the two systems. 

5. For a suitable variation of the initial density distribu- 
tion E ~ r~ d of the disk, with d between 1.5 and 2 the 
spectral energy distributions of DQ Tau and AK Sco 
could be fitted well to the observational data. 

Future models of this kind need to include a more detailed 
structure of the stars as well as a more detailed model of 
the accretion process. The numerical resolution we have 
reached by now is so fine that the individual stars in close 
binaries are already resolved by usually about 100 grid- 
cells within low-res calculations. 

To obtain a detailed comparison with observed light 
curves in particular the phase dependence a more elabo- 
rate model of the radiative loss of the very thin material 
in the circumstellar disks is required. In that case 3d sim- 
ulations may be necessary which are beyond present day 
computational resources. 
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